*COMDECK COM1
      COMMON/C1/Q(500,6,60)
      COMMON/C2/X(500,60),Y(500,60),Z(500,60)
      COMMON /C3/ ALP,RE,FSMACH,GAMMA,JMAX,KMAX,LMAX,JWMAX,LWF,ND,LE,FN
     1,NOPT,JT,ETAE,ITMAX,LMAXBC,L1BC,SMU,DT,PR,TW,KREADY,IRGRID,ICRUDE,
     2  ITOT,INWALL,LW,KW,JREADX,IWRIT,JLW,JKW,KLW,FSP,FST,JL1WL,JL1WU,
     3   KMAXBC,JK1WU,JK1WL
      COMMON/C4/IGEOM,XMAX,XZERO,YMAX,X1,X2,X3,X4,Z1,Z2,Z3,Z4,ZT,RT
      COMMON /C5/EMC(60),UC(60),PC(60),TC(60),UE(60),DELT(60),DELTZ(60)
      COMMON /C6/RI(60),PTI(60),HTI(60),VROUI(60),PTOT(500),HTOT(500),
     1    VOU(600),WOU(600),RL(60)
      COMMON/C7/YKR(60),TH(60)
      COMMON/C8/XX(62,3)
*DECK CRUDIC
      SUBROUTINE CRUDIC
*CALL COM1
      EMC2(B)=1.+FLOAT(IJ)*B**.5
      FUNC(A)=4.*EMC2(A)*(1.+(GAMMA-1.)/2.*EMC2(A))
      IF(NOPT.EQ.3) GO TO 60
C
C       NOPT=2
C
      LW1=LW+1
      DO 55 K=1,KMAX
C     SET UP INITIAL CONDITIONS AT INFLOW PLANE J=1
C
C     INSIDE NOZZLE FOR NOPT=2 WE WILL USE INITQ RESULTS
C
C      OUTSIDE THE NOZZLE
C
      DO 35 L=LW1,LMAX
      KL=(L-1)*ND+K
      PCTOT=(FSP/PTOT(KL))**((GAMMA-1.)/GAMMA)
      VLS=2.*HTOT(KL)*(1.-PCTOT)/(GAMMA-1.)
      VL=SQRT(VLS)
      Q(KL,1,1)=PTOT(KL)/HTOT(KL)*(FSP/PTOT(KL))**(1.0/GAMMA)
      Q(KL,2,1)=Q(KL,1,1)*VL/SQRT(1.+VOU(KL)**2+WOU(KL)**2)
      Q(KL,3,1)=Q(KL,2,1)*VOU(KL)
      Q(KL,4,1)=Q(KL,2,1)*WOU(KL)
      T=HTOT(KL)*PCTOT
      Q(KL,5,1)=Q(KL,1,1)*(T/(GAMMA*(GAMMA-1.))+VLS/2.)
   35 CONTINUE
C
C     SET UP INITIAL CONDITIONS IN PLANES J > 0
C
C     OUTSIDE THE NOZZLE L > LW  (REGION 1 AND 2)
C
      DO 50 L=LW1,LMAX
      KL=(L-1)*ND+K
      DO 50 J=2,JMAX
      IF (L.EQ.LMAX) GO TO 45
      DO 40 N=1,5
   40 Q(KL,N,J)=Q(KL,N,1)
      GO TO 50
C
C     FREESTREAM BOUNDARY L=LMAX
C
   45 Q(KL,1,J)=FSP/FST
      Q(KL,2,J)=Q(KL,1,J)*(FSMACH*SQRT(FST))
      Q(KL,3,J)=0.
      Q(KL,4,J)=0.
      Q(KL,5,J)=FSP*(1./(GAMMA*(GAMMA-1.))+.5*FSMACH**2)
   50 CONTINUE
   55 CONTINUE
C     ENSURE VELOCITIES =0 AT WALLS
      DO 57 J=1,JMAX
      DO 57 L=LW,LW1
      DO 57 K=1,KMAX
      KL=(L-1)*ND+K
      Q(KL,2,J)=0.
      Q(KL,3,J)=0.
      Q(KL,4,J)=0.
      IF(L.EQ.LW) GO TO 57
      Q(KL,1,J)=PTOT(KL)/HTOT(KL)
      Q(KL,5,J)=PTOT(KL)/(GAMMA*(GAMMA-1.))
   57 CONTINUE
      GO TO 999
   60 CONTINUE
C
C       NOPT=3
C
      JKLW=MIN0(JKW,JLW)
      IF(KW.EQ.0) KMW=KMAX
      IF(KW.NE.0) KMW=KW
      IF(LW.EQ.0) LMW=LMAX
      IF(LW.NE.0) LMW=LW
C
C  CALCULATE LOCAL STATIC PRESSURE
C
      T0=1.0
      P0=1.0
      DO 65 I=1,2
      IF(JT.EQ.JMAX.AND.I.EQ.1) GO TO 65
      IF(I.EQ.1) IJ=+1
      IF(I.EQ.2) IJ=-1
      J=JT
      PHIJ=0.
      XXX=X(1,J)
      CALL GEOM(XXX,AJ,ZZT)
      EMC(J)=1.
   63 TC(J)=T0/(1.+(GAMMA-1.)/2.*EMC(J)**2)
      PC(J)=P0*TC(J)**(GAMMA/(GAMMA-1.))
      IF(J.EQ.JMAX) GO TO 65
      IF(J.EQ.1) GO TO 65
      J=J+IJ
      AB=AJ
      XXX=X(1,J)
      CALL GEOM(XXX,AJ,ZZT)
      DELPSI=ALOG(AJ/AB)
      PHIHAT=PHIJ+DELPSI*FUNC(PHIJ)
      PHIJ=.5*(PHIJ+PHIHAT+DELPSI*FUNC(PHIHAT))
      EMC(J)=EMC2(PHIJ)**.5
      GO TO 63
   65 CONTINUE
      DO 100 J=1,JMAX
      DO 70 K=1,KMAX
      DO 70 L=1,LMAX
      KL=(L-1)*ND+K
      P=FSP
      IF(K.LE.KW.AND.L.LE.LW.AND.J.LT.JKLW) P=PC(J)
      CALL QCALC(J,K,L,HTOT(KL),PTOT(KL),P)
   70 CONTINUE
C
C     ZERO VELOCITY AT ALL WALLS
C
      LWP1= LW + 1
      KWP1= KW + 1
      IF(J.LT.JKLW) GO TO 80
      IF(J.GT.JKW) GO TO 75
      DO 72 L=1,LWP1
      DO 72 K=KW,KWP1
      KL=(L-1)*ND+K
      P=FSP
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   72 CONTINUE
   75 IF(J.GT.JLW) GO TO 90
      DO 77 K=1,KWP1
      DO 77 L=LW,LWP1
      KL=(L-1)*ND+K
      P=FSP
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   77 CONTINUE
      GO TO 90
   80 IF(KW.EQ.0) GO TO 86
      DO 85 K=KW,KWP1
      IF(K.EQ.KW) P=PC(J)
      IF(K.EQ.KW+1) P=FSP
      IF(K.EQ.KW) LS=LW
      IF(K.EQ.KW+1) LS=LW+1
      IF(LW.EQ.0) LS=LMAX
      DO 85 L=1,LS
      KL=(L-1)*ND+K
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   85 CONTINUE
   86 IF(LW.EQ.0) GO TO 90
      DO 87 L=LW,LWP1
      IF(L.EQ.LW) P=PC(J)
      IF(L.EQ.LW+1) P=FSP
      IF(L.EQ.LW) KS=KW
      IF(L.EQ.LW+1) KS=KW+1
      IF(KW.EQ.0) KS=KMAX
      DO 87 K=1,KS
      KL=(L-1)*ND+K
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   87 CONTINUE
   90 IF(J.LT.JL1WL) GO TO 95
      IF(J.GT.JL1WU) GO TO 95
C
C   WEDGE PLUG AT L=1
C
      DO 92 K=1,KMW
      L=1
      IF(J.GE.JKLW) P=FSP
      IF(J.LT.JKLW) P=PC(J)
      KL=(L-1)*ND+K
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   92 CONTINUE
   95 CONTINUE
      IF(JK1WU.EQ.0) GO TO 99
      IF(J.LT.JK1WL) GO TO 99
      IF(J.GT.JK1WU) GO TO 99
C
C    WEDGE PLUG AT K=1
C
      DO 97 L=1,LMW
      K=1
      IF(J.LT.JKLW) P=PC(J)
      IF(J.GE.JKLW) P=FSP
      KL=(L-1)*ND+K
      CALL QZERO(J,K,L,HTOT(KL),PTOT(KL),P)
   97 CONTINUE
   99 CONTINUE
C
  100 CONTINUE
  999 RETURN
      END
*DECK DELI
      FUNCTION DELI(X)
      DELI=.025*RWFUNC(X)
      RETURN
      END
*DECK DELOUT
      FUNCTION DELOUT(X)
      DELOUT=.01*(RMAX(X)-RWOUT(X))
      RETURN
      END
*DECK DKMET
      SUBROUTINE DKMET(J,K,L,KL,XK,YK,ZK)
*CALL COM1
C     INTERMEDIATE WALL IN K DIRECTION
      DY2 = .5
      KP = KL+1
      KR = KL-1
      IF(KW.LE.0.OR.K.LT.KW.OR.K.GT.KW+1.OR.J.GT.JKW) GO TO 200
      IF(LW.GT.0.AND.L.GT.LW) GO TO 200
      IF(K.EQ.KW) GO TO 700
      IF(K.EQ.KW+1) GO TO 500
  200 IF(K.NE.1) GO TO 600
C     TEST FOR PLANAR SYMMETRY
      IF(KPLANE.EQ.0) GO TO 250
      IF(KMAX.EQ.1) GO TO 1000
C     TEST FOR INTERMEDIATE WALL NORMAL TO K=1 SURFACE
  250 IF(LW) 300,300,450
  300 IF(JK1WL.LE.J.AND.J.LE.JK1WU) GO TO 500
C     SYMMETRY
  400 XK=0.0
      ZK=0.0
      YK=2.0*(Y(KP,J)-Y(KL,J))*DY2
      RETURN
C     TEST FOR WALL NORMAL TO K=1 SURFACE
  450 IF(L-LW) 300,300,400
C     FOWARD DIFFERECNE
  500 FAC=2.0
      KR=KL
      GO TO 900
  600 IF(K.NE.KMAX) GO TO 800
C     TEST FOR PLANAR SYMMETRY
      IF(KPLANE.EQ.0) GO TO 650
      IF(KMAX.EQ.1) GO TO 1000
  650 IF(KMAXBC.LT.3.OR.KMAXBC.GT.4) GO TO 700
C     SYMMETRY
      XK=0.0
      YK=0.0
      ZK=0.0
      KR=KR
      IF(KMAXBC.EQ.3) YK=2.0*(Y(KL,J)-Y(KR,J))*DY2
      IF(KMAXBC.EQ.4) ZK=2.0*(Z(KL,J)-Z(KR,J))*DY2
      RETURN
C     BACKWARD DIFFERENCE
  700 KP=KL
      FAC=2.0
      GO TO 900
C     CENTRAL DIFFERENCE
  800 FAC=1.0
  900 KR=KR
      KP=KP
      XK=(X(KP,J)-X(KR,J))*DY2*FAC
      YK=(Y(KP,J)-Y(KR,J))*DY2*FAC
      ZK=(Z(KP,J)-Z(KR,J))*DY2*FAC
      RETURN
C     PLANAR SYMMETRY
 1000 XK=0.0
      YK=1.0
      ZK=0.0
      RETURN
      END
*DECK DLMET
      SUBROUTINE DLMET(J,K,L,KL,XL,YL,ZL)
*CALL COM1
      COMMON/AXISYM/LAXIS
C     INTERMEDIATE WALL IN L DIRECTION
      DZ2 = .5
      LP = KL+ND
      LR = KL-ND
      IF(LW.LE.0.OR.L.LT.LW.OR.L.GT.LW+1.OR.J.GT.JLW) GO TO 100
      IF(KW.GT.0.AND.K.GT.KW) GO TO 100
      IF(L.EQ.LW) GO TO 700
      IF(L.EQ.LW+1) GO TO 500
  100 IF(L.NE.1) GO TO 600
C     AXIS OF SYMMETRY
      IF(LAXIS.NE.1) GO TO 150
      XL=0.0
      YL=Y(LP,J)-Y(KL,J)
      ZL=Z(LP,J)-Z(KL,J)
      GO TO 999
  150 CONTINUE
C     TEST FOR INTERMEDIATE WALL NORMAL TO L=1 SURFACE
      IF(KW) 200,200,400
C     TEST FOR WALL AT L=1
  200 IF(JL1WL.LE.J.AND.J.LE.JL1WU) GO TO 500
C     SYMMETRY
  300 XL=0.0
      YL=0.0
      ZL=2.0*(Z(LP,J)-Z(KL,J))*DZ2
      GO TO 999
C     TEST FOR WALL AT L=1
  400 IF(K-KW) 200,200,300
C     FORWARD DIFFERENCE
  500 LR=KL
      FAC=2.0
      GO TO 900
  600 IF(L.NE.LMAX) GO TO 800
      IF(LMAXBC.LT.3.OR.LMAXBC.GT.4) GO TO 700
C     SYMMETRY
      XL=0.0
      YL=0.0
      ZL=0.0
      LR=LR
      IF(LMAXBC.EQ.4) ZL=2.0*(Z(KL,J)-Z(LR,J))*DZ2
      IF(LMAXBC.EQ.3) YL=2.0*(Y(KL,J)-Y(LR,J))*DZ2
      GO TO 999
C     BACKWARD DIFFERENCE
  700 LP=KL
      FAC=2.0
      GO TO 900
C     CENTRAL DIFFERENCE
  800 FAC=1.0
  900 LR=LR
      LP=LP
      XL=(X(LP,J)-X(LR,J))*DZ2*FAC
      YL=(Y(LP,J)-Y(LR,J))*DZ2*FAC
      ZL=(Z(LP,J)-Z(LR,J))*DZ2*FAC
  999 RETURN
      END
*DECK GEOM
      SUBROUTINE GEOM(XIN,A,ZIN)
*CALL COM1
      DIMENSION AR(32),B(32)
      DATA PI2/1.57079632/
      DATA IFIRST/0/
      ZFUNC(XXX)=ZT+XXX*XXX/(RT+(RT*RT-XXX*XXX)**.5)
      IF(NOPT.EQ.2) GO TO 120
      IF(NOPT.EQ.3) GO TO 150
      IF(IGEOM.EQ.2) GO TO 35
C     2D COSINE NOZZLE
      XT=0.
      ZI=2.0
      ZE=1.5
      ZT=1.0
C     NOZZLE HEIGHT Z AND AREA A AT STATION X RELATIVE TO THROAT XT
      IF(XIN-XT)10,10,20
  10  ZIN=ZT+.5*(ZI-ZT)*(1.+COS(PI2*(XIN-XZERO)))
      GO TO 30
  20  ZIN=ZT+.5*(ZE-ZT)*(1.+COS(PI2*(XMAX-XIN)))
C     2D NOZZLE AREA
  30  A=ZIN
      GO TO 999
   35 CONTINUE
      IF(IFIRST.NE.0) GO TO 38
      ZT=.5388
      RT=1.0777
      X1=-2.275
      X2=-.4095
      X3=.0228
      X4=2.275
      Z1=1.3859
      Z2=ZFUNC(X2)
      Z3=ZFUNC(X3)
      Z4=.5868
      IFIRST=1
   38 CONTINUE
      IF(XIN-X1) 40,40,50
   40 ZIN=Z1
      GO TO 100
   50 IF(XIN-X2) 60,60,70
   60 ZIN=Z1+((Z2-Z1)/(X2-X1))*(XIN-X1)
      GO TO 100
   70 IF(XIN-X3) 80,80,90
   80 ZIN=ZFUNC(XIN)
      GO TO 100
   90 ZIN=Z3+((Z4-Z3)/(X4-X3))*(XIN-X3)
  100 A=ZIN*YMAX
      GO TO 999
C
C     NOPT=2
C
  120 CONTINUE
      A=PI2/2.*RWFUNC(XIN)**2
      GO TO 999
C
C      NOPT=3
C
  150 KT=KMAX
      IF(KW.GT.0) KT=KW
      LT=LMAX
      IF(LW.GT.0) LT=LW
      DO 153 J=1,JMAX
      IF(X(1,J).NE.XIN) GO TO 153
      JP=J
      GO TO 154
  153 CONTINUE
  154 CONTINUE
      DO 160 L=1,LT
      KLO=(L-1)*ND
      DO 155 K=1,KT
      KL=KLO+K
      M=K
      CALL XXM(M,L,JP,JP)
      DA=SQRT(XX(JP,1)*XX(JP,1)+XX(JP,2)*XX(JP,2)+XX(JP,3)*XX(JP,3))
  155 AR(K)=DA
      CALL QDRTR(B(L),AR,1.,1,KT)
  160 CONTINUE
      CALL QDRTR(A,B,1.,1,LT)
  999 RETURN
      END
*DECK INFLOW
      SUBROUTINE INFLOW
*CALL COM1
      ITOT=LMAX
      INWALL =LW
      XXX=X(1,1)
      RWI=RWFUNC(XXX)*ZT
      RWO=RWOUT(XXX)*ZT
      GAMI=GAMMA/(GAMMA-1.)
      DELI=.05
      DELO=1.0
      DO 35 L=1,LW
      RI(L)=RL(L)
      RR=(RWI-RI(L))/DELI
      IF(RR.GE.1) GO TO 20
      IF(RR.GT.3.E-3) GO TO 15
      UUE=145.3*RR
      GO TO 30
   15 UUE=RR**(1./7.)
      GO TO 30
   20 UUE=1.
   30 CONTINUE
      HTI(L)=1.
      VROUI(L)=0.
      PTI(L)=(1.+(GAMMA-1)/2.*EMC(1)**2*(1.-UUE*UUE))**(-GAMI)
      WRITE(6,100) L,RI(L),PTI(L),HTI(L),VROUI(L)
  100 FORMAT(22H L,RI,PTI,HTI,VROUI= ,I3,1P4E12.5)
   35 CONTINUE
      LW1=LW+1
      DO 70 L=LW1,LMAX
      RI(L)=RL(L)
      VROUI(L)=0.
      UUO=(RI(L)-RWO)/DELO
      IF(UUO.GT.1.) GO TO 50
      UUO=UUO**(1./7.)
      GO TO 60
   50 UUO=1.
   60 SAVE=(GAMMA-1.)/2.*FSMACH**2
      HTI(L)=FST*(1.+SAVE)
      PTI(L)=FSP*((1.+SAVE)/(1.+SAVE*(1-UUO*UUO)))**GAMI
      WRITE(6,100) L,RI(L),PTI(L),HTI(L),VROUI(L)
   70 CONTINUE
      RETURN
      END
*DECK INITQ
      SUBROUTINE INITQ
C     INITIALIZE FLOW VARIABLES Q(KL,N,J) AT GRID POINTS
*CALL COM1
      DIMENSION ZRL(32),PT(32),HT(32),VROU(32)
      DIMENSION ETA(13),FP(13),G(13)
      DIMENSION TE(32)
      DATA ETA/0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0/
      DATA FP/0.0,.16586,.32979,.48652,.62977,.75073,.84605,.91255,
     1                               .95552,.97928,.99155,.99682,.99898/
      DATA G/0.0,.03978,.16422,.35812,.60951,.87940,1.14133,1.3557,
     1                          1.51632,1.61644,1.67446,1.70189,1.71424/
      EMC2(B)=1.+FLOAT(IJ)*B**.5
      FUNC(A)=4.*EMC2(A)*(1.+(GAMMA-1.)/2.*EMC2(A))
C     INITIAL DATA FOR FLAT PLATE BOUNDARY LAYER(NOPT=0)
      IF (NOPT.GT.0) GO TO 40
      DO 10 I=1,13
   10 ETA(I)=ETA(I)/6.0
      DO 35 L=LWF,LMAX
      ZETA=(L-1.0)/(LMAX-1.0)
      ZETAE=(LE-1.)/(LMAX-1.0)
      ZFAC=(FN+1.)*(FN**(ZETA/(1.-ZETAE))-1.)/(FN**(1./(1.-ZETAE))-1.)
      IF(ZFAC.EQ.0.) EZFAC=1.
      IF(ZFAC.EQ.0.) GO TO 12
      EZFAC=EXP(-(6.1*ZFAC)**2)
   12 G4=1.71424*FSMACH*EZFAC/2.0
      G5=1.0/(GAMMA*(GAMMA-1.0))
      DO 35 K=1,KMAX
      KL=(L-1)*ND+K
      IE=2
      DO 35 J=1,JWMAX
      ETAJ=Z(KL,J)*SQRT(RE/X(KL,J))/ETAE
      IF (ETAJ.GT.1.0) GO TO 30
      DO 15 I=IE,13
      IF (ETA(I).GT.ETAJ) GO TO 20
   15 CONTINUE
      I=12
   20 IF(I.EQ.13) I=12
      IE=I
      CALL TRPOL8(ETAJ,ETA,FP,I-1,FPJ)
      CALL TRPOL8(ETAJ,ETA,G,I-1,GJ)
      U=FSMACH*FPJ
      V=0.0
      W=FSMACH*GJ/(2.0*SQRT(X(KL,J)*RE))
      Q(KL,1,J)=1.0/(1.0+.5*(GAMMA-1.0)*FSMACH**2*(1.0-(U/FSMACH)**2))
      IF (TW.LE.0.0) GO TO 25
      TB=TW+(U/FSMACH)*(1.0-TW+(GAMMA-1.0)/2.0*FSMACH**2*(1.0-U/FSMACH))
      Q(KL,1,J)=1.0/TB
   25 Q(KL,2,J)=Q(KL,1,J)*U
      Q(KL,3,J)=Q(KL,1,J)*V
      Q(KL,4,J)=Q(KL,1,J)*W
      Q(KL,5,J)=G5+.5*(Q(KL,2,J)**2+Q(KL,3,J)**2+Q(KL,4,J)**2)/Q(KL,1,J)
      GO TO 35
   30 Q(KL,1,J)=1.0
      Q(KL,2,J)=FSMACH
      Q(KL,3,J)=0.0
      Q(KL,4,J)=G4/SQRT(X(KL,J)*RE)
      IF(FSMACH.LE.1.0) Q(KL,4,J)=1.71424*FSMACH/(2.0*SQRT(X(KL,J)*RE))
      Q(KL,5,J)=G5+.5*(Q(KL,2,J)**2+Q(KL,3,J)**2+Q(KL,4,J)**2)/Q(KL,1,J)
   35 Q(KL,6,J)=1.0
      GO TO 999
   40 CONTINUE
      IF (NOPT.EQ.3) GO TO 80
      LMW=LMAX
      IF(LW.GT.0) LMW=LW
      LE=(LMAX+1)/2
      IF(LW.GT.0) LE=(LW+1)/2
C     INITIAL DATA FOR INTERIOR OF 2D NOZZLE
C     SET STAGNATION PRESSURE AND TEMPERATURE
      T0=1.0
      P0=1.0
C     1-D ISENTROPIC CHOKED FLOW
      DO 75 I=1,2
      IF (JT.EQ.JMAX.AND.I.EQ.1) GO TO 75
      IF(I.EQ.1) IJ=+1
      IF(I.EQ.2) IJ=-1
      J=JT
      PHIJ=0.
      XXX=X(1,J)
      CALL GEOM(XXX,AJ,ZTT)
      EMC(J)=1.
C     AVERAGE FLOW VARIABLES
   45 TC(J)=T0/(1.+(GAMMA-1.)/2.*EMC(J)**2)
      PC(J)=P0*TC(J)**(GAMMA/(GAMMA-1.))
      UC(J)=EMC(J)*TC(J)**.5
      DO 70 L=1,LMW
      DO 70 K=1,KMAX
      KL=(L-1)*ND+K
      IF (L.GT.LE) GO TO 50
      UL=UC(J)
      TL=TC(J)
      GO TO 60
   50 UL=UC(J)*FLOAT(LMW-L)/FLOAT(LMW-LE)
      IF (TW.GT.0) GO TO 55
      TL=TC(J)*(1.+(GAMMA-1.)/2.*EMC(J)**2*(1.-(UL/UC(J))**2))
      GO TO 60
   55 TL=TW+(UL/UC(J))*(-TW+TC(J)*(1.+(GAMMA-1.)/2.*EMC(J)**2*(1.-UL/
     1  UC(J))))
   60 VL=0.
      CALL XYZXI(KL,J,XJ,YJ,ZJ)
      WL=UL*ZJ/XJ
      IF (NOPT.LT.2) GO TO 65
      IF(K.EQ.1) VR=UL*ZJ/XJ
      VL=VR*SIN(TH(K))
      WL=VR*COS(TH(K))
   65 RHOL=PC(J)/TL
      EL=PC(J)/(GAMMA*(GAMMA-1.))+RHOL*((UL*UL+VL*VL+WL*WL)/2.)
      Q(KL,1,J)=RHOL
      Q(KL,2,J)=RHOL*UL
      Q(KL,3,J)=RHOL*VL
      Q(KL,4,J)=RHOL*WL
      Q(KL,5,J)=EL
      Q(KL,6,J)=1.0
   70 CONTINUE
      IF (J.EQ.JMAX) GO TO 75
      IF (J.EQ.1) GO TO 75
      J=J+IJ
      AB=AJ
      XXX=X(1,J)
      CALL GEOM(XXX,AJ,ZTT)
      DELPSI=ALOG(AJ/AB)
      PHIHAT=PHIJ+DELPSI*FUNC(PHIJ)
      PHIJ=.5*(PHIJ+PHIHAT+DELPSI*FUNC(PHIHAT))
      EMC(J)=EMC2(PHIJ)**.5
      GO TO 45
   75 CONTINUE
      IF(NOPT.EQ.2.AND.ITOT.EQ.0.AND.IGEOM.EQ.0) CALL INFLOW
   80 CONTINUE
      IF (ITOT.EQ.0) GO TO 140
      DO 95 I=1,ITOT
      IF (LW.GT.0) GO TO 85
      RI(I)=RI(I)/RI(ITOT)
      GO TO 95
   85 IF (I.GT.INWALL) GO TO 90
      RI(I)=RI(I)/RI(INWALL)
      GO TO 95
   90 IF(I.EQ.INWALL+1) RISAVE=RI(INWALL+1)
      RI(I)=(RI(I)-RISAVE)/(RI(ITOT)-RISAVE)
   95 CONTINUE
      DO 115 L=1,LMAX
      KL=(L-1)*ND+1
      LL=LMAX
      IF (LW.EQ.0) GO TO 105
      LL=LW
      IF (L.LE.LW) GO TO 105
      KLW1=LW*ND+1
      KLMAX=(LMAX-1)*ND+1
      IF (NOPT.GT.1) GO TO 100
      ZRL(L)=(Z(KL,1)-Z(KLW1,1))/(Z(KLMAX,1)-Z(KLW1,1))
      GO TO 115
  100 ZRL(L)=(RL(L)-RL(LW+1))/(RL(LMAX)-RL(LW+1))
      GO TO 115
  105 IF (NOPT.GT.1) GO TO 110
      KL2=(LL-1)*ND+1
      ZRL(L)=Z(KL,1)/Z(KL2,1)
      GO TO 115
  110 ZRL(L)=RL(L)/RL(LL)
  115 CONTINUE
C
C     INTERPOLATE RADIAL DISTRIBUTION OF PRESSURE AND TOTAL ENTHALPY
C     TO ACTUAL GRID
C
      DO 130 L=1,LMAX
      IF(L.GT.LW) GO TO 131
      DO 125 I=1,INWALL
      C1=ZRL(L)-RI(I)
      C2=RI(I+1)-ZRL(L)
      IF (C1*C2) 125,120,120
  120 PT(L)=(C1*PTI(I+1)+C2*PTI(I))/(C1+C2)
      HT(L)=(C1*HTI(I+1)+C2*HTI(I))/(C1+C2)
      VROU(L)=(C1*VROUI(I+1)+C2*VROUI(I))/(C1+C2)
      GO TO 130
  125 CONTINUE
  131 INW1=INWALL+1
      DO 134 I=INW1,ITOT
      C1=ZRL(L)-RI(I)
      C2=RI(I+1)-ZRL(L)
      IF (C1*C2) 134,132,132
  132 PT(L)=(C1*PTI(I+1)+C2*PTI(I))/(C1+C2)
      HT(L)=(C1*HTI(I+1)+C2*HTI(I))/(C1+C2)
      VROU(L)=(C1*VROUI(I+1)+C2*VROUI(I))/(C1+C2)
      GO TO 130
  134 CONTINUE
  130 CONTINUE
      DO 135 L=1,LMAX
      DO 135 K=1,KMAX
      KL=(L-1)*ND+K
      PTOT(KL)=PT(L)
      HTOT(KL)=HT(L)
      VOU(KL)=VROU(L)*SIN(TH(K))
      WOU(KL)=VROU(L)*COS(TH(K))
  135 CONTINUE
  140 CONTINUE
      IF (ICRUDE.NE.0) GO TO 180
      IF (ITOT.EQ.0) GO TO 150
      DO 145 K=1,KMAX
      DO 145 L=1,LMAX
      KL=(L-1)*ND+K
      PCTOT=(PC(1)/PTOT(KL))**((GAMMA-1.)/GAMMA)
      VLS=2.*HTOT(KL)*(1.-PCTOT)/(GAMMA-1.)
      VL=SQRT(VLS)
      Q(KL,1,1)=PTOT(KL)/HTOT(KL)*(PC(1)/PTOT(KL))**(1.0/GAMMA)
      Q(KL,2,1)=Q(KL,1,1)*VL/SQRT(1.+VOU(KL)**2+WOU(KL)**2)
      Q(KL,3,1)=Q(KL,2,1)*VOU(KL)
      Q(KL,4,1)=Q(KL,2,1)*WOU(KL)
      T=HTOT(KL)*PCTOT
      Q(KL,5,1)=Q(KL,1,1)*(T/(GAMMA*(GAMMA-1.))+VLS/2.)
  145 CONTINUE
  150 CONTINUE
      IF(KREADY.GT.0) CALL SIDWIC
C     BOUNDARY LAYER THICKNESS FOR ADIABATIC B.C.
      WRITE (6,155)
  155 FORMAT(//T4,1HJ,T11,3HUCJ,T24,3HTCJ,T37,3HUEJ,T50,3HTEJ,
     1 T62,5HDELTJ,T74,6HDELTZJ)
C
      KL=(LMAX-1)*ND+1
      SJ=0.
      XXX=X(1,1)
      CALL GEOM(XXX,AJ,ZTT)
      DELT(1)=FN/(1.+FN)*ZT
      DO 175 J=1,JMAX
      CALL XYZXI(KL,J,XJ,YJ,ZJ)
      UE(J)=UC(J)*(1.+(ZJ/XJ)**2)**.5
      TE(J)=1.-(GAMMA-1.)/2.*UE(J)**2
      FJ=(2.+TE(J))/(UE(J)**5*TE(J)**((2.-GAMMA)/(GAMMA-1.)))
      IF(J.EQ.1)F1=FJ
      GJ=UE(J)**9*TE(J)**((4.-3.*GAMMA)/(GAMMA-1.))*(XJ*XJ+ZJ*ZJ)**.5
      IF (J.EQ.1) GO TO 160
      SJ=SJ+.5*(GJ+GJ1)
      DELT(J)=FJ*((DELT(1)/F1)**2+(4./(3.*RE))*SJ)**.5
  160 DELTZ(J)=DELT(J)*(1.+(ZJ/XJ)**2)**.5
      GJ1=GJ
      WRITE (6,170)J,UC(J),TC(J),UE(J),TE(J),DELT(J),DELTZ(J)
  170 FORMAT(2X,I2,6(2X,1PE11.4))
  175 CONTINUE
      GO TO 185
  180 CALL CRUDIC
  185 CONTINUE
  999 RETURN
      END
*DECK INPUT
      SUBROUTINE INPUT
*CALL COM1
      DATA LWF/1/
C     NOPT
C            0   EXTERNAL 2-D BOUNDARY LAYERS
C            1   INTERNAL FLOW IN 2-D NOZZLES
C            2   AXIZYMMETRIC NOZZLES  AXIS IS AT Y=Z=0  (THE CARTESIAN
C                    L=1,K=1,KMAX
C                PURE INTERNAL FLOW OR COMBINED EXTERNAL/INTERNAL FLOW
C                FOR KW>0 OR LW>0  SET ICRUDE=1
C                FOR JL1WU > 0 SET ICRUDE=1
C            3   ARBITRARY NONAXISYMMETRIC NOZZLES
C                FOR KW > 0 OR LW > 0 SET ICRUDE =1
C
C     IRGRID
C            0   COMPUTE THE GRID
C            1   READ IN THE GRID FROM A FILE PREPARED BY RGRID
C
C     ICRUDE
C            0   USE THE CODE INITQ TO GET THE INITIAL Q'S
C            1   USE THE CODE CRUDIC TO GET THE INITIAL Q'S
      READ (5,10)NMAX,JMAX,KMAX,LMAX,LAMIN,INVISC,J1BC,JMAXBC,KPLANE,K1B
     1C,JK1WL,JK1WU,KW,JKW,KMAXBC
      READ (5,10)L1BC,JL1WL,JL1WU,LW,JLW,LMAXBC,NRST,IWRIT,NGRI,NP,KVIS
     1,LVIS,KLVIS,INFLT,ISUTH,NROUT
      READ (5,15)DT,FSMACH,RMACH,RE,PR,RTDEGK,FSP,FST
      READ (5,15)GAMMA,RMUEXP,TW,CNBR,DTFAC,RM,SMU,OMEGA
   10 FORMAT(16I5)
   15 FORMAT(8F10.0)
      READ (5,10)NOPT,JREADX,IGEOM,KREADY,IRGRID,ICRUDE,ITOT,INWALL
      READ (5,15)XZERO,XMAX,FN,YMAX
      IF (JREADX.GT.0) READ (5,15)(X(1,J),J=1,JREADX)
      IF (KREADY.GT.0) READ (5,15)(YKR(K),K=1,KREADY)
      IF (ITOT.EQ.0) GO TO 25
      DO 20 I=1,ITOT
      READ (5,15)RI(I),PTI(I),HTI(I),VROUI(I)
   20 CONTINUE
   25 ALP=0.
      ITMAX=NMAX
      WRITE (6,30)NMAX,JMAX,KMAX,LMAX,LAMIN,INVISC,J1BC,JMAXBC,KPLANE,K1
     1BC,JK1WL,JK1WU,KW,JKW,KMAXBC
      WRITE (6,35)L1BC,JL1WL,JL1WU,LW,JLW,LMAXBC,NRST,IWRIT,NGRI,NP,KVIS
     1,LVIS,KLVIS,INFLT,ISUTH,NROUT
      WRITE (6,40)DT,FSMACH,RMACH,RE,PR,RTDEGK,FSP,FST
      WRITE (6,45)GAMMA,RMUEXP,TW,CNBR,DTFAC,RM,SMU,OMEGA
   30 FORMAT(122H1    NMAX    JMAX    KMAX    LMAX    LAMIN  INVISC   J1
     1BC   JMAXBC  KPLANE   K1BC    JK1WL   JK1WU    KW      JKW   KMAXB
     1C     /16I8)
   35 FORMAT(//132H     L1BC    JL1WL   JL1WU    LW      JLW   LMAXBC
     1 NRST   IWRIT   NGRI     NP      KVIS    LVIS   KLVIS   INFLT   IS
     2UTH    NROUT         /16I8)
   40 FORMAT(//7X,2HDT,10X,6HFSMACH,10X,5HRMACH
     1,12X,2HRE,13X,2HPR,11X,6HRTDEGK,10X,3HFSP,13X,3HFST,/1P8E15.7)
   45 FORMAT(//5X,5HGAMMA,10X,6HRMUEXP,11X,2HTW,12X,4HCNBR,10X,5HDTFAC,
     113X,2HRM,12X,3HSMU,11X,5HOMEGA/1P8E15.7///)
      WRITE (6,50)NOPT,JREADX,IGEOM,KREADY,IRGRID,ICRUDE,ITOT,INWALL
   50 FORMAT(/  66H     NOPT  JREADX   IGEOM  KREADY  IRGRID  ICRUDE
     1ITOT  INWALL              /16I8    )
      WRITE (6,55)XZERO,XMAX,FN,YMAX
   55 FORMAT(/5X,5HXZERO,10X,4HXMAX,13X,2HFN,12X,4HYMAX,
     1     /1P4E15.7)
      IF (JREADX.GT.0) WRITE (6,60)(X(1,J),J=1,JREADX)
   60 FORMAT(/8H  X(J)=  (8F10.5))
      IF (KREADY.GT.0) WRITE (6,65)(YKR(K),K=1,KREADY)
   65 FORMAT(8H  Y(K)=   (8F10.5))
      IF (ITOT.EQ.0) GO TO 85
      WRITE (6,70)
   70 FORMAT(T4,1HI,T16,2HRI,T28,3HPTI,T41,3HHTI,T53,5HVROUI)
      DO 80 I=1,ITOT
      WRITE (6,75)I,RI(I),PTI(I),HTI(I),VROUI(I)
   75 FORMAT(2X,I2,6X,4(2X,1PE11.4))
   80 CONTINUE
   85 CONTINUE
      IF(IRGRID.EQ.0) GO TO 110
      IF(NOPT.NE.3) GO TO 110
      WRITE(6,89)
   89 FORMAT(T8,1HL,T15,1HK,T24,4HPTOT,T38,4HHTOT,T52,3HVOU,T65,3HWOU)
      NKL=KMAX*LMAX
      DO 100 I=1,NKL
      READ(5,90) L,K,PTOT((L-1)*KMAX+K),HTOT((L-1)*KMAX+K),
     1  VOU((L-1)*KMAX+K),WOU((L-1)*KMAX+K)
   90 FORMAT(2I5,4F10.0)
      WRITE(6,95) L,K,PTOT((L-1)*KMAX+K),HTOT((L-1)*KMAX+K),
     1  VOU((L-1)*KMAX+K),WOU((L-1)*KMAX+K)
   95 FORMAT(4X,2(I5,2X),4(2X,1PE11.4))
  100 CONTINUE
  110 CONTINUE
      RETURN
      END
*DECK LAGRAN
      SUBROUTINE LAGRAN (XX,X,FT)
      DIMENSION X(4),FT(4),A(4),B(4),C(4),D(4),E(4)
      DO 5 I=1,4
      A(I)=XX-X(I)
      B(I)=X(1)-X(I)
      C(I)=X(2)-X(I)
      D(I)=X(3)-X(I)
    5 E(I)=X(4)-X(I)
      IF (XX.GT.X(2)) GO TO 7
      FT(4)=0.
      FT(1)=(A(2)*A(3))/(B(2)*B(3))
      FT(2)=(A(1)*A(3))/(C(1)*C(3))
      FT(3)=(A(1)*A(2))/(D(1)*D(2))
      RETURN
    7 IF (XX.LT.X(3)) GO TO 8
      FT(1)=0.
      FT(2)=(A(3)*A(4))/(C(3)*C(4))
      FT(3)=(A(2)*A(4))/(D(2)*D(4))
      FT(4)=(A(2)*A(3))/(E(2)*E(3))
      RETURN
    8 FT(1)=(A(2)*A(3)*A(4))/(B(2)*B(3)*B(4))
      FT(2)=(A(1)*A(3)*A(4))/(C(1)*C(3)*C(4))
      FT(3)=(A(1)*A(2)*A(4))/(D(1)*D(2)*D(4))
      FT(4)=(A(1)*A(2)*A(3))/(E(1)*E(2)*E(3))
      RETURN
      END
*DECK MAIN
      PROGRAM NOZLIC(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,TAPE2=ICFILE,
     1  TAPE3=GPLUG1)
C     MAIN ROUTINE TO SET UP GRID AND INITIALIZE FLOW VARIABLES
      CALL INPUT
      CALL NDLPTS
      CALL INITQ
      CALL OUTPUT
      STOP
      END
*DECK NDLPTS
      SUBROUTINE NDLPTS
*CALL COM1
      DATA LWF/1/
      PI=4.*ATAN(1.)
      ND2=KMAX*LMAX
      ZT=0.
      IF(NOPT.GT.0.AND.FN.NE.0) FN=(1.0-FN)/FN
      ETAE=6.0*(1.0+0.07*FSMACH*(0.2+FSMACH))
      JWMAX=JMAX
      ND=KMAX
      IF (IRGRID.EQ.0) GO TO 20
      DO 15 J=1,JMAX
      READ(3) X(1,J),(Y(KL,J),Z(KL,J),KL=1,ND2)
      DO 10 N=2,ND2
   10 X(N,J)=X(1,J)
   15 CONTINUE
      GO TO 85
   20 CONTINUE
      DO 80 K=1,KMAX
      YJKL=0.
      TH(1)=0.
      IF (KMAX.EQ.1) GO TO 25
      TH(K)=PI/2.*FLOAT(K-1)/FLOAT(KMAX-1)
      YJKL=YMAX*FLOAT(K-1)/FLOAT(KMAX-1)
   25 CONTINUE
      DO 80 L=LWF,LMAX
      KL=(L-1)*ND+K
      DO 80 J=1,JWMAX
      IF (JREADX.GT.0) GO TO 35
      X(KL,J)=XZERO+FLOAT(J-1)/FLOAT(JMAX-1)*(XMAX-XZERO)
      GO TO 40
   35 IF (KL.EQ.1) GO TO 40
      X(KL,J)=X(1,J)
   40 CONTINUE
      CALL STRCH(J,K,L,ZFAC)
      IF (NOPT.LT.2) GO TO 55
      XXX=X(KL,J)
      IF (LW.GT.0.AND.L.GT.LW) GO TO 45
      RJL=RWFUNC(XXX)*ZFAC
      GO TO 50
   45 RJL=RWOUT(XXX)+(RMAX(XXX)-RWOUT(XXX))*ZFAC
   50 Y(KL,J)=RJL*SIN(TH(K))
      Z(KL,J)=RJL*COS(TH(K))
          IF(K.EQ.1 .AND. J.EQ.1) RL(L)=RJL
      GO TO 80
   55 IF (KREADY.GT.0) GO TO 60
      Y(KL,J)=YJKL
      GO TO 65
   60 IF(K.LE.KREADY) Y(KL,J)=YKR(K)
   65 IF (NOPT.NE.0) GO TO 70
      ZZ=ETAE*SQRT(X(KL,J)/RE)
      GO TO 75
   70 XXX=X(KL,J)
      CALL GEOM(XXX,AA,ZZ)
   75 Z(KL,J)=ZZ*ZFAC
   80 CONTINUE
      IF(KREADY.GT.0 .AND. KREADY.LT.KMAX) CALL YSTRCH(YKR(KREADY))
   85 CONTINUE
C   SEARCH FOR THROAT
      DO 90 J=1,JMAX
      IF (X(1,J).NE.0) GO TO 90
      JT=J
      GO TO 95
   90 CONTINUE
      JT=1
   95 CONTINUE
      WRITE(6,100) JT
  100 FORMAT(2X,21HNOZZLE THROAT AT JT=  I3)
      IF(NOPT.GE.1) CALL SCALE
      RETURN
      END
*DECK OUTPUT
      SUBROUTINE OUTPUT
*CALL COM1
      IT=0
      GD=GAMMA*(GAMMA-1.)
      TAU=0.0
      CNBR=0.0
      NK=0
      ND2=KMAX*LMAX
      WRITE(2) KMAX,JMAX,LMAX,ITMAX,LMAXBC,L1BC,FSMACH,GAMMA,RE,SMU,DT
     1                                                      ,ALP,CNBR,PR
      WRITE(2) IT,TAU,DT,NK
      WRITE(2) ((X(KL,J),KL=1,ND2),J=1,JMAX),
     1         ((Y(KL,J),KL=1,ND2),J=1,JMAX),
     2         ((Z(KL,J),KL=1,ND2),J=1,JMAX)
      WRITE(2) (((Q(KL,N,J),KL=1,ND2),N=1,6),J=1,JMAX)
      IF (IWRIT.EQ.0) GO TO 35
      DO 30 J=1,JMAX
      DO 30 K=1,KMAX
      WRITE (6,15)J,K
   15 FORMAT(1H0,2X,2HJ=,I3,2X,2HK=,I3,2X,1HL,6X,1HX,11X,1HY,11X,1HZ
     1  ,6X,6HR/RREF,5X,6HU/AREF,5X,6HV/AREF,5X,6HW/AREF,5X,6HT/TREF,
     1 5X,6HP/PREF,5X,3HENT)
      DO 25 L=1,LMAX
      KL = (L-1)*ND+K
      R = Q(KL,1,J)
      RR = 1./R
      U = Q(KL,2,J)*RR
      V = Q(KL,3,J)*RR
      W = Q(KL,4,J)*RR
      E = Q(KL,5,J)
      S2=0.0
      IF(ABS(U).GT.1.0E-17) S2=S2+U**2
      IF(ABS(V).GT.1.0E-17) S2=S2+V**2
      IF(ABS(W).GT.1.0E-17) S2=S2+W**2
      PP = GD*(E-.5*R*S2)
      TT=PP*RR
      ENT = PP/(ABS(R))**GAMMA
      WRITE (6,20)L,X(KL,J),Y(KL,J),Z(KL,J),R,U,V,W,TT,PP,ENT
   20 FORMAT(1H ,14X,I3,10F11.6)
   25 CONTINUE
   30 CONTINUE
   35 CONTINUE
C     IF(NOPT.EQ.1) WRITE(3)EMC,UC,PC,TC,UE,DELT,DELTZ
      RETURN
      END
*DECK QCALC
      SUBROUTINE QCALC(J,K,L,HT,PT,P)
*CALL COM1
      KL=(L-1)*ND+K
      CALL XYZXI(KL,J,XJ,YJ,ZJ)
      T=HT*(P/PT)**((GAMMA-1.)/GAMMA)
      RHO=P/T
      VLS=2./(GAMMA-1.)*ABS(HT-T)
      V=SQRT(VLS)
      R=SQRT(XJ*XJ+YJ*YJ+ZJ*ZJ)
      Q(KL,1,J)=RHO
      Q(KL,2,J)=RHO*V*XJ/R
      Q(KL,3,J)=RHO*V*YJ/R
      Q(KL,4,J)=RHO*V*ZJ/R
      Q(KL,5,J)=RHO*(T/(GAMMA*(GAMMA-1.))+VLS/2.)
      RETURN
      END
*DECK QDRTR
      SUBROUTINE QDRTR(ENTGRL,ENTGRD,DLT,IL,IU)
      DIMENSION ENTGRD(32)
C   1-D TRAPEZOIDAL QUADRATURE. COMPUTES DEFINITE INTEGRAL,ENTGRL,
C   OF INTEGRAND,ENTGRD(I),BETWEEN LIMITS IL,IU.
C
      ENTGRL=0.5*(ENTGRD(IL)+ENTGRD(IU))
       IF((IU-IL).LT.2) RETURN
      ILP=IL+1
      IUM=IU-1
      DO 1 I=ILP,IUM
    1 ENTGRL=ENTGRL+ENTGRD(I)
      ENTGRL=ENTGRL*DLT
      RETURN
      END
*DECK QZERO
      SUBROUTINE QZERO(J,K,L,HT,PT,P)
*CALL COM1
      KL=(L-1)*ND+K
      T=HT*(P/PT)**((GAMMA-1.)/GAMMA)
      RHO=P/T
      Q(KL,1,J)=RHO
      Q(KL,2,J)=0.
      Q(KL,3,J)=0.
      Q(KL,4,J)=0.
      Q(KL,5,J)=RHO*T/(GAMMA*(GAMMA-1.))
      RETURN
      END
*DECK RMAX
      FUNCTION RMAX(XIN)
*CALL COM1
      DIMENSION XINPUT(28),RIN(28)
      DATA XINPUT/-34.5,-29.5,-24.5,-19.5,-15.5,-12.5,-10.1,-8.4,-7.0,
     1   -5.842,-4.953,-4.064,-3.048,
     2  -2.032,-1.016,0.,1.08,2.17,3.26,4.35,5.44,6.53,7.62,8.7,
     3  10.,11.5,13.,14.5/
      DATA RIN/28*12.62/
      DATA IFIRST/0/
      IF(ZT.EQ.0.) GO TO 5
      IF(IFIRST.NE.0) GO TO 5
      IFIRST=1
      DO 2 J=1,JMAX
      XINPUT(J)=XINPUT(J)/ZT
      RIN(J)=RIN(J)/ZT
    2 CONTINUE
    5 CONTINUE
      DO 10 J=1,JMAX
      IF(XIN.NE.XINPUT(J)) GO TO 10
      JFIND=J
      GO TO 20
   10 CONTINUE
      GO TO 990
   20 CONTINUE
      RMAX=RIN(JFIND)
      RETURN
  990 WRITE(6,999)
  999 FORMAT(2X,35H***ERROR IN RMAX FUNCTION ROUTINE**)
      STOP
      END
*DECK RWFUNC
      FUNCTION RWFUNC(XIN)
*CALL COM1
      DIMENSION XINPUT(28),RIN(28)
      DATA XINPUT/-34.5,-29.5,-24.5,-19.5,-15.5,-12.5,-10.1,-8.4,-7.0,
     1   -5.842,-4.953,-4.064,-3.048,
     2  -2.032,-1.016,0.,1.08,2.17,3.26,4.35,5.44,6.53,7.62,8.7,
     3  10.,11.5,13.,14.5/
      DATA RIN/4*6.4315,6.336,6.194,5.996,5.820,5.653,5.436,5.019,4.585,
     1  4.229,3.993,3.856,13*3.81/
      DATA IFIRST/0/
      IF(ZT.EQ.0.) GO TO 5
      IF(IFIRST.NE.0) GO TO 5
      IFIRST=1
      DO 2 J=1,JMAX
      XINPUT(J)=XINPUT(J)/ZT
      RIN(J)=RIN(J)/ZT
    2 CONTINUE
    5 CONTINUE
      DO 10 J=1,JMAX
      IF(XIN.NE.XINPUT(J)) GO TO 10
      JFIND=J
      GO TO 20
   10 CONTINUE
      GO TO 990
   20 CONTINUE
      RWFUNC=RIN(JFIND)
      RETURN
  990 WRITE(6,999) XIN,XINPUT
  999 FORMAT(2X,37H***ERROR IN RWFUNC FUNCTION ROUTINE**/(8E12.5))
      STOP
      END
*DECK RWOUT
      FUNCTION RWOUT(XIN)
*CALL COM1
      DIMENSION XINPUT(28),RIN(28)
      DATA XINPUT/-34.5,-29.5,-24.5,-19.5,-15.5,-12.5,-10.1,-8.4,-7.0,
     1   -5.842,-4.953,-4.064,-3.048,
     2  -2.032,-1.016,0.,1.08,2.17,3.26,4.35,5.44,6.53,7.62,8.7,
     3  10.,11.5,13.,14.5/
      DATA RIN/4*7.62,7.54624,7.38481,7.1898,
     1  7.01605,6.85065,6.69854,6.57232,
     2  6.43787,6.27408,6.09943,5.91386,5.71732,
     3  5.49628,5.26045,5.01173,4.75002,4.47521,
     4   4.18718,6*3.812/
      DATA IFIRST/0/
      IF(ZT.EQ.0) GO TO 5
      IF(IFIRST.NE.0) GO TO 5
      IFIRST=1
      DO 2 J=1,JMAX
      XINPUT(J)=XINPUT(J)/ZT
      RIN(J)=RIN(J)/ZT
    2 CONTINUE
    5 CONTINUE
      DO 10 J=1,JMAX
      IF(XIN.NE.XINPUT(J)) GO TO 10
      JFIND=J
      GO TO 20
   10 CONTINUE
      GO TO 990
   20 CONTINUE
      RWOUT=RIN(JFIND)
      RETURN
  990 WRITE(6,999)
  999 FORMAT(2X,36H***ERROR IN RWOUT FUNCTION ROUTINE**)
      STOP
      END
*DECK SCALE
      SUBROUTINE SCALE
*CALL COM1
      DIMENSION R(13)
      EQUIVALENCE (XMAX,R(1))
      LL=LMAX
      IF(LW.GT.0) LL=LW
      KL=(LL-1)*ND+1
      ZT=Z(KL,JT)-Z(1,JT)
      ZTT=ZT
      IF(NOPT.NE.1) IOP=2
      IF(NOPT.EQ.1) IOP=13
      DO 10 J=1,JMAX
      DO 10 L=1,LMAX
      DO 10 K=1,KMAX
      KL=(L-1)*ND+K
      X(KL,J)=X(KL,J)/ZT
      Y(KL,J)=Y(KL,J)/ZT
      Z(KL,J)=Z(KL,J)/ZT
   10 CONTINUE
      DO 20 I=1,IOP
      R(I)=R(I)/ZTT
   20 CONTINUE
      RETURN
      END
*DECK SIDWIC
      SUBROUTINE SIDWIC
C
C      SIDE WALL INITIAL CONDITIONS
*CALL COM1
      KE=(KMAX+1)/2
      DO 70 J=1,JMAX
      DO 65 L=1,LMAX
      KLE=(L-1)*ND+KE
      U0=Q(KLE,2,J)/Q(KLE,1,J)
      DO 65 K=KE,KMAX
      UK=U0*(FLOAT(KMAX-K)/FLOAT(KMAX-KE))
      KL=(L-1)*ND+K
      IF (TW.GT.0) GO TO 55
      TL=TC(J)*(1.+(GAMMA-1.)/2.*EMC(J)**2*(1.-(UK/UC(J))**2))
      GO TO 60
   55 TL=TW+(UK/UC(J))*(-TW+TC(J)*(1.+(GAMMA-1.)/2.*EMC(J)**2*(1.-UK/
     1  UC(J))))
   60 CONTINUE
      RHOL=PC(J)/TL
      EL=PC(J)/(GAMMA*(GAMMA-1.))+RHOL*UK*UK/2.
      Q(KL,1,J)=RHOL
      Q(KL,2,J)=RHOL*UK
      Q(KL,3,J)=0.
      Q(KL,4,J)=0.
      Q(KL,5,J)=EL
   65 CONTINUE
   70 CONTINUE
      RETURN
      END
*DECK STRCH
      SUBROUTINE STRCH(J,K,L,ZFAC)
*CALL COM1
      DATA IFIRST/0/
      IF(IFIRST.NE.0) GO TO 5
      IFIRST=1
      IFN=0
      IF(FN.NE.0.) IFN=1
    5 KL=(L-1)*KMAX+K
      IF(IFN.NE.0) GO TO 12
      XXX=X(KL,J)
      IF(LW.NE.0 .AND. L.GT.LW) GO TO 11
      FN=DELI(XXX)/RWFUNC(XXX)
      GO TO 12
   11 FN=DELOUT(XXX)/(RMAX(XXX)-RWOUT(XXX))
   12 IF(LW.GT.0) GO TO 14
      LE=(LMAX+1)/2
      ZETA=FLOAT(L-1)/FLOAT(LMAX-1)
      ZETAE=FLOAT(LE-1)/FLOAT(LMAX-1)
      GO TO 30
   14 IF(L.GT.LW) GO TO 20
      LE=(LW+1)/2
      ZETA=FLOAT(L-1)/FLOAT(LW-1)
      ZETAE=FLOAT(LE-1)/FLOAT(LW-1)
      GO TO 30
   20 LE=(LMAX-LW+1)/2+LW+1
      ZETA=1.-FLOAT(L-LW-1)/FLOAT(LMAX-LW-1)
      ZETAE=FLOAT(LE-LW-1)/FLOAT(LMAX-LW-1)
   30 ZFAC=(FN+1.)*(FN**(ZETA/(1.-ZETAE))-1.)/(FN**(1./(1.-ZETAE))-1.)
      IF(NOPT.LT.1) GO TO 99
      ZFAC=ZFAC/(FN+1.)
      IF(L.GT.LW.AND.LW.GT.0) ZFAC=1.-ZFAC
   99 RETURN
      END
*DECK TRPOL8
      SUBROUTINE TRPOL8(XX,X,Y,I,YY)
      DIMENSION X(1),A(4),Y(1)
      IF(I.EQ.1) GO TO 100
      CALL LAGRAN(XX,X(I-1),A)
      YY=A(1)*Y(I-1)+A(2)*Y(I)+A(3)*Y(I+1)+A(4)*Y(I+2)
      RETURN
  100 YY=Y(1)+(Y(2)-Y(1))*(XX-X(1))/(X(2)-X(1))
      RETURN
      END
*DECK XXM
      SUBROUTINE XXM(M,L,J1,J2)
*CALL COM1
C
C  XI METRICS FORMED FOR A K,L LINE IN J
C
C
C  SYMMETRY
C
      K = M
      KL = (L-1)*ND+K
      DO 10 J = J1,J2
      CALL DKMET(J,K,L,KL,XK,YK,ZK)
      CALL DLMET(J,K,L,KL,XL,YL,ZL)
      XX(J,1) = YK*ZL-ZK*YL
      XX(J,2) = ZK*XL-XK*ZL
      XX(J,3) = XK*YL-YK*XL
   10 CONTINUE
      RETURN
      END
*DECK XYZXI
      SUBROUTINE XYZXI(KL,J,XJ,YJ,ZJ)
*CALL COM1
      DX2=.5
      JP=J+1
      JR=J-1
C     XI DERIVATIVES OF X,Y,Z
      IF(J.EQ.1) GO TO 50
      IF(J.EQ.JMAX) GO TO 51
      XJ = (X(KL,JP)-X(KL,JR))*DX2
      YJ = (Y(KL,JP)-Y(KL,JR))*DX2
      ZJ = (Z(KL,JP)-Z(KL,JR))*DX2
      GO TO 70
   50 J1 = J+1
      XJ=X(KL,J1)-X(KL,J)
      YJ=Y(KL,J1)-Y(KL,J)
      ZJ=Z(KL,J1)-Z(KL,J)
      GO TO 70
   51 J1 = J-1
      XJ=X(KL,J)-X(KL,J1)
      YJ=Y(KL,J)-Y(KL,J1)
      ZJ=Z(KL,J)-Z(KL,J1)
   70 CONTINUE
      RETURN
      END
*DECK YSTRCH
      SUBROUTINE YSTRCH(YKREAD)
C      ROUTINE TO EXPONENTIALLY STRETCH THE Y GRID FROM Y(KREADY+1) TO Y
C
*CALL COM1
      FD(W)=((W-1.)*EXP(W)+1.)/(W*W)
      F(W)=(EXP(W)-1.)/W-1./BETA
      DELE=1./(KMAX-KREADY)
      DO 35 J=1,JMAX
      KLMAX=(LMAX-1)*ND+1
      KLM=(LMAX-2)*ND+1
      DELZ=Z(KLMAX,J)-Z(KLM,J)
      BETA=DELZ/((YMAX-YKREAD)*DELE)
C
C     NEWTON ITERATION FOR Y STRETCHING COEFFICIENT
      W1=20.
      N=0
   10 W2=W1-F(W1)/FD(W1)
      N=N+1
      ERR=ABS((W2-W1)/W2)
      IF (ERR.LT.1.E-5) GO TO 20
      W1SAVE=W1
      W1=W2
      IF (N.LT.20) GO TO 10
C
C     FAILURE TO CONVERGE
C
      WRITE (6,15)W1SAVE,W2,ERR,J,N,YMAX,YKREAD,DELE,DELZ,BETA
   15 FORMAT(81H****FAILURE TO CONVERGE IN NEWTON ITERATION FOR Y STRETC
     1HING COEFFICIENT OMEGA***   /2X,3HW1=E10.3,4H W2=E10.3,5H ERR=,
     2  E10.3,3H J=I2,3H N=I2,/2X,6H YMAX= E10.3,8H YKREAD=E10.3,
     3  6H DELE=  E10.3,  6H DELZ=E10.3,6H BETA=E10.3)
      GO TO 40
   20 CONTINUE
      OMEGA=W2
      OMEINV=1./(1.-EXP(-OMEGA))
      KR1=KREADY+1
        DO 30 K=KR1,KMAX
        ETAK=FLOAT(K-KREADY)/FLOAT(KMAX-KREADY)
        EXPOET=(1.-EXP(-OMEGA*ETAK))*OMEINV
          DO 25 L=1,LMAX
          KL=(L-1)*ND+K
          Y(KL,J)=YKREAD+(YMAX-YKREAD)*EXPOET
   25     CONTINUE
   30   CONTINUE
   35 CONTINUE
      RETURN
C*** ERROR STOP
   40 STOP
      END
